4  Using Meta-regressions to Compare CDD across Species or Sites

4.1 Overview

Following from Section 4. How does CDD vary across species, abiotic gradients or in time? in the main text, here, we demonstrate one approach to comparing the strength of CDD across species using a meta-analysis framework. The same approach can be used to compare CDD across sites, plots, or any other unit of interest, as long as it is possible to generate reliable estimates of the strength of CDD (suitable sample sizes, etc.). As noted in the main text, “Correct propagation of uncertainty in CDD estimates requires meta-regressions in frequentist or Bayesian frameworks. Weighted regressions (e.g., lm, lmer, gam) can also estimate how CDD varies e.g., across latitude or between species with different life-history strategies, but incorrectly estimate the associated uncertainty.

We use a subset of the BCI seedling data (Comita et al. 2023) of only 30 species to compare how species abundance is related to the strength of CDD at the species level. This is the same dataset used in section 2. This analysis is for demonstration purposes only, and biological conclusions should not be made about the results, given this is only a small subset of the data.

An advantage of using meta-regressions over simple weighted regressions is that the models are able to simultaneously account for uncertainty in species-specific CDD estimates as well as systematic differences in species’ CDD when regressed against a predictor via the inclusion of random effects. A simple weighted regression, on the other hand, assumes that there is only one error for which the relative strength is known when regressing the estimates against a predictor. The latter can lead to incorrect weighting of species in the metaregression.

In this tutorial, we use relative average marginal effect (𝑟𝐴𝑀𝐸) calculated in the previous chapter as our response variable, calculated separately for each species. 𝑟𝐴𝑀𝐸 in this case estimates the relative increase in the probability of annual mortality with the addition of one new conspecific neighbor, while keeping total densities constant. Positive numbers indicate a relative increase in mortality with an increase in conspecific density, a signature of NCDD. In principle, any metric of the strength of CDD can be used, though care must be taken to ensure that the metrics are comparable across species and sites (see main text for more information).

We use the popular metafor package to fit the meta-regression models.

Note

The following code is adapted from the latitudinalCNDD repository by Lisa Hülsmann.

4.2 Load libraries and data

Code
# Load libraries
library(tidyr)
library(dplyr)
library(readr)
library(ggplot2)
library(here)
library(metafor)

# Load in species abundances for BCI data subset
abundances <- read_csv(
here("./data/BCI seedling data - 30 species - abundance 2023_05_18.csv")
  )

# Load marginal effects calculations from previous section
load(here("./data/mortality.Rdata"))

# Subset down to just equilibrium change
rAME <- rAME %>%
  filter(change == "equilibrium")

# Join marginal effects and abundance data
rAME <- left_join(rAME, abundances, by = c("sp" = "spp"))

# Add in average abundance for 'rare' species
rAME$abundance[rAME$sp == "data_deficient_seedling"] <- 
  mean(abundances$abundance[abundances$spp %in% 
                              nsp$spp[nsp$data_deficient]])

# Log transform abundance to use in models
rAME$log_abundance <- log(rAME$abundance)

Let’s take a quick look at the data set we’ll be working with:

Code
head(rAME, n = 10)
       term    estimate   std.error offset change.value      change     sp
1  con_dens  0.12190548 0.122096554      1          + 1 equilibrium ACALDI
2  con_dens  0.46332386 0.140120748      1          + 1 equilibrium AEGICE
3  con_dens  0.02026319 0.003624488      1          + 1 equilibrium BEILPE
4  con_dens  0.42931069 0.087277501      1          + 1 equilibrium CAPPFR
5  con_dens  0.04121409 0.042571758      1          + 1 equilibrium CECRIN
6  con_dens  0.10440629 0.241547627      1          + 1 equilibrium CORDLA
7  con_dens  0.11658356 0.200294935      1          + 1 equilibrium DAVINI
8  con_dens  0.01113154 0.123539647      1          + 1 equilibrium DESMAX
9  con_dens  0.23375216 0.227252222      1          + 1 equilibrium HEISCO
10 con_dens -0.06498043 0.061086207      1          + 1 equilibrium JUSTGR
   abundance log_abundance
1        169      5.129899
2        153      5.030438
3       5574      8.625868
4       1604      7.380256
5         51      3.931826
6        248      5.513429
7         27      3.295837
8         55      4.007333
9         56      4.025352
10       136      4.912655

4.3 Reformat data for model fitting

First, we use the ‘escalc’ function in the metafor package to essentially repackage our data frame into a format used in the meta-regression model fitting. Since we already calculated our effect size (𝑟𝐴𝑀𝐸), we just pass through the 𝑟𝐴𝑀𝐸 estimate and corresponding standard error using the ‘GEN’ option for the ‘measure’ argument, rather than calculating an effect size within the ‘escalc’ function.

Code
# Reformat data for model fitting
# Set measure to generic, which passes the observed effect sizes or 
# outcomes via the yi argument and the corresponding sampling 
# variances via the vi argument (or the standard errors via the sei 
# argument) to the function.
    dat_meta = metafor::escalc(measure = "GEN", 
                               yi = estimate, # observed outcomes
                               sei = std.error, # standard errors
                               slab = sp, # label for species
                               data = rAME)

4.4 Fit meta-regression model

Next, we use the ‘rma’ function to fit a meta-regression model, where 𝑟𝐴𝑀𝐸 is our response variable (renamed as yi in the previous step) and log species abundance is our predictor. While not shown here, it is possible to fit mixed effects meta-regression models with the ‘rma.mv’ function. We suggest consulting the extensive documentation for the metafor package for further details.

Code
# Fit model
metamod = metafor::rma(yi = yi,
                       vi = vi,
                       mods = ~ log_abundance,
                       method = "REML",
                       data = dat_meta)

4.6 Plot the estimated of rAME for all species with a forest plot

In this case, a forest plot shows the estimates of the strength of CDD for individual species, here ordered by least to most abundant going from top to bottom.

Code
forest(metamod, 
       header = "Species", 
       xlab = "rAME",
       order = log_abundance)

4.7 Model diagnostics

The plot method displays model diagnostics of the meta-regression model (more info here) in addition to a forest plot.

Code
plot(metamod)